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Sparse matter is characterized by regions with low electron density and its understanding calls for 
methods to accurately calculate both the van der Waals (vdW) interactions and other bonding. Here 
we present a first-principles density functional theory (DFT) study of a layered oxide (V2O5) bulk 
structure which shows charge voids in between the layers and we highlight the role of the vdW forces 
in building up material cohesion. The result of previous first-principles studies involving semilocal 
approximations to the exchange-correlation functional in DFT gave results in good agreement with 
experiments for the two in-plane lattice parameters of the unit cell but overestimated the parameter 
for the stacking direction. To recover the third parameter we include the nonlocal (dispersive) 
vdW interactions through the vdW-DF method [Dion et al., Phys. Rev. Lett. 92, 246401 (2004)] 
testing also various choices of exchange fiavors. We find that the transferable first-principle vdW- 
DF calculations stabilizes the bulk structure. The vdW-DF method gives results in fairly good 
agreement with experiments for all three lattice parameters. 

PACS numbers: 31.15.E-,71.15.Mb,71.15.Nc,61.50.Lt 



I. INTRODUCTION 

Vanadium is one of the most abundant metals on earth 
and it is found in about 150 different minerals. Catalysts 
based on vanadium oxides are widely used in the produc- 
tion of chemicals and in the reduction of environmental 
pollution, in particular the vanadium pentoxide (V2O5) 
form is extensively used.l^In recent years V2O5 has also 
been used with intercalating Li ions for high-capacity 
solid-state batteries.'^JEl Xhe applied-research and indus- 
trial focus on catalysis and batteries involving V2O5 has 
led to a substantial amount of atomic-scale theory studies 
focusing on the V2O5 bulkP^^and its surfaces. I^ti^ 

The bulk of V2O5 has a stacked structure. Previ- 
ous atomic-scale calculations that were based on den- 
sity functional theory (DFT) with the popular semilo- 
cal generalized gradient approximations (GGA) have of- 
ten failed in arriving at the experimentally known value 
of the lattice parameter in the stacking direction of the 
layers.l^'l^ This is believed to happen because GGA does 
not describe the van der Waals (vdW) forces,^^^ and be- 
cause these are expected to play an important role in 
binding the layers in V2O5 bulk.E^In several GGA-based 
V2O5 surface or vacancy studies this deficiency of GGA 
was either ignored or worked around by imposing the ex- 
perimentally obtained lattice parameter or unit cell vol- 
ume. Other groups have employed the semi-empirical 
method DFT-D (ReL \TQ for adding the vdW forces.li^l 

In this paper we use the first -princ iples vdW density- 
functional approach, vdW-DF,'2SEil to determine the 
V2O5 bulk structure. We find that the vdW forces are in- 
deed important for binding the structure, and we analyse 
the binding within and mainly between the layers. We 
here focus solely on the common a-V205 bulk structure 
and ignore the more exotic 7-V2O5 structure.^ 

The vdW-DF method has a different philosophy and 
aim and some advantages over the DFT-D method. The 



DFT-D approach to include the vdW forces is an atom- 
centered pair-potentials method. The vdW-DF method 
is directly based on the electron response. The vdW-DF 
method correctly describes the interaction as arising in 
the tails of the electron distribution, not at the atomic 
centers. Unlike DFT-D, vdW-DF provides a framework 
which is well suited to include effects of image planes 
The DFT-D omission of image planes can result in incon- 
sistencies of the description across a range of distances.!^ 
This effect could be important in materials like V2O5 
where the corrugation is large, and where as a result 
both relatively small and larger distances contribute to 
the vdW interaction simultaneously. 

The outline is as follows. In Sec. II we describe the 
structure of V2O5 as known from experiments. In Sec. 
Ill we describe the computational methods used, both 
the self-consistent (sc) GGA calculations that provide 
the electron density, and the postprocess procedure that 
takes this density as input for obtaining the vdW contri- 
bution. Sec. IV is devoted to the discussion of our results 
and Sec. V includes a summary. 



II. MATERIAL STRUCTURE 

The (q;-)V2 05 bulk has orthorhombic symmetry and a 
layered structure: each vanadium atom is connected to 
five oxygen atoms to create pyramids that share their cor- 
ners in building a double chain. The chains are connected 
along the edges to form layers that are then stacked to 
form the bulk structure (Figure [T]). 

There are three structurally different oxygen atoms in 
each layer. One is coordinated to one vanadium atom, 
the second is found in a bridging position and the third 
has a threefold coordinated position. The binding of the 
atoms inside a layer is strong whereas the interactions 
that keep the layers stacked are weaker, resulting in V2O5 
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FIG. 1: The structure of V2O5 bulk (q phase). The box shows 
the primitive unit cell that contains two formula units. The 
14 atoms of the unit cell are shown, as well as a number of the 
atoms of the neighboring unit cells. V atoms are represented 
by large spheres and O atoms by small spheres. The labels a, 
b, and c indicate the vectors that span the orthorhombic unit 
cell. The cleavage plane is perpendicular to the c direction. 
Figure created using xcrysden;-'* 



being easily cleaved. 



III. COMPUTATIONAL METHODS 

For determining the bulk structure of V2O5 we use 
the first-principles nonlocal functional vdW-DF within 
DFT.I^ We calculate the vdW-DF energies in a post- 
GGA procedure.l^Sl 

The main interest here is on determining the optimal 
inter-layer spacing in V2O5 (which is the same as the 
unit cell side length c, Figure [T]) and the binding en- 
ergy, defined as the energy gained by stacking layers of 
V2O5. We calculate the vdW-DF total energy £;vdW-DF 
for a number of values of c and find the minimum of 
the energy curve. For each point on the energy curve 
the procedure requires first a sc-GGA calculation, from 
which the GGA-based total energy E^f^ and the sc- 
GGA electron density n are evaluated. Then n is used for 
evaluating the long-range correlation contribution arising 
from the vdW interactions, Following the system- 

atic proced ure desc ribed in more detail in several other 
publicationSj'^^IzSIlT) g^j^^ summarized below, we combine 
the sc-GGA and nonlocal results to obtain ^;vdW-DF 

This section contains three parts explaining the com- 
putational methods used. The first part deals with the 
sc-GGA calculations used as a basis for obtaining n and 
some terms of the energy i<;'<'dW-DF rpj^^ second part con- 
tains a short summary of the scheme used for the calcula- 
tion of Ef and the vdW-DF total energy £;vdW-DF^ 
the third part briefly explains the effect and importance 
of choice of exchange functional for use in vdW-DF. 



A. The sc-GGA calculations 



The sc-GGA calculations are carried out using the 
plane-wave code—- dacapo with ultrasoft pseudopoten- 
tials (USPP). We describe V2O5 by an orthorhombic unit 
cell containing two formula units, periodically repeated 
in all directions. The energy cutoff for the expansion of 
the wave functions is set to 500 eV. The sc-GGA cal- 
culations are carried out using the PBE (Ref. 29) flavor 
of the GGA for the exchange and correlation functional. 
The fast Fourier transform (FFT) grid is chosen such 
as to have a distance less than 0.12 A between nearest- 
neighbor grid points. The choice of this relatively dense 
electron density grid is important for the quality of the 
subsequent evaluation of the nonlocal correlation contri- 
bution. The Brillouin zone of the unit cell is sampled 
according to the Monkhorst-Pack scheme by means of a 
2x4x4 fc-point sampling. 

Our results are converged with respect to the wave 
function cut off, the number of k points needed, and a 
number of code-internal parameters. By lowering the 
number of k points in the c direction we notice an in- 
significant 0.001 eV) change in the minimum value of 
the sc-GGA total energy curve. 

Since the intraplanar bonds have an ionic and partly 
covalent nature^ we assume the layers to remain rigid 
when changing the interlayer spacing c. This means that 
the atoms are kept fixed in their positions (found at the 
PBE binding distance) relative to the layer as the dis- 
tance between the layers is changed, with fixed atom- 
atom distances within the layers. By analysing the resid- 
ual Hellmann-Feynman forces on the atoms, calculated 
from the GGA-based electron density, we find that this 
is a reasonable approximation. 

In this paper we define the binding energy Eb of V2O5 
to be the energy gained by moving together isolated lay- 
ers of V2O5 to form the bulk structure. Each unit cell in 
the bulk contains two formula units in one layer, periodi- 
cally repeated, and the reference ( "layers far apart" ) cal- 
culation is carried out with vacuum added on both sides 
of the layer. The reference calculation is for a unit cell 
that has the c direction side length four times that of the 
original unit cell, periodically repeated in all directions.'^ 
We use the same reference unit cell for the nonlocal cor- 
relation (S"') reference calculations and this imposes ad- 
ditional constraints on the construction of the reference 
unit cell, as described further below. 

To supplement our analysis of the electron density with 
a better description of the core electrons we carry out a 
set of additional GGA calculations using the all-electron 
DFT code GPAwEHthat is based on projector augmented 
waveJ^ (PAW). In these calculations we use settings of 
computational parameters as close as possible to the pa- 
rameters we use in DACAPO. 
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B. vdW density functional 

We use the vdW-DF scheme to include the vdW inter- 
actions in a systematic manner. The correlation energy 
Ec is splilP^ into a nearly-local part E^ and a part that 
includes the most nonlocal interactions i?"', 

E,^E",+Ef. (1) 

The splitting of the correlation contributions makes it 
possible to employ different approximations for each 
term. In a homogeneous system the term E^ is the cor- 
relation E^"^^ obtained from the local density approxi- 
mation (LDA), and in generaP^l we approximate E^ by 
^LDA rpj^g j^ni vanishes for a homogeneous system. It 
describes the coupling through the electrodynamic field, 
the dispersion interaction. The difference in contri- 
butions provides a description of the interaction which 
acts across large distances. The difference is not much 
influenced by local variations in the electron density. 
Rather, it is susceptible to the more coarse-grained re- 
sponse of the environment. The form of is derived in 
Ref.EOland is 

Ef[n]^^J j drdr' n{v)<j){Yy)n{r') (2) 

where (j) \s a. kernel, explicitly stated in Ref. [20l 

The electron density n that enters ^ is taken from the 
sc-GGA calculations and is described on a grid (the FFT 
grid of the sc-GGA calculation) . Our code for evaluating 
i?"' is not periodic. However, as described in detail in 
Refs. I22I31I351 the natural periodicity within the bulk 
unit cell is easily represented by explicitly including a 
number of the periodic images of the unit cell, with the 
nearby grid points carefully described and the grid points 
far away only described via a coarse version of the grid. 

Similar to the energy contributions from the sc-GGA 
calculations the term i?"' must be evaluated relative to 
a reference calculation with the V2O5 layers "far apart". 
As described above, this is here achieved with a unit cell 
four times the original unit cell in the c direction, with 
the V2O5 layer placed close to the middle of the unit 
cell. For this reference calculation of E^^ no periodicity 
in the c direction is assumed. The layer must be placed 
such that locally, with respect to the (FFT) grid points 
on which n is described, the nuclei maintain the same 
positions in the bulk and the reference calculation. This 
means that the (FFT) grid must be chosen such that the 
bulk and the reference calculations have the exact same 
grid point separation; in this study we need precisely 
four times the number of FFT points in the reference 
calculation compared to the original bulk calculation. 

The combination of correlation terms in ([T]) avoids dou- 
ble counting of correlation contributions. Using the cor- 
relation energy E^. from ([T]) the total energy in the vdW- 
DF scheme is 

^vdW-DF _ ^GGA _ ^GGA _^ ^LDA _|_ ^nl 



where E^^^ is the total energy from the sc-GGA calcula- 
tion and E^'^^ and E^^^ the GGA respective LDA cor- 
relation calculated from the sc-GGA electron density n. 
The layer binding energy E^ (per unit cell) we define as 
the energy gained by moving the layers in V2O5 together 
accordion-like,'^ 

TP ^pvdW-DF TpvdW-DF\ ( ,\ 

Eb ~ -{E^nWa -^ref )■ W 

C. Exchange functionals 

With every choice of approximation for correlation 
functional comes a need for choosing a suitable flavor of 
exchange functional. In previous work we have chosen to 
combine the correlation Ec of ([T]) with the exchange from 
the GGA flavor revPBE.^- This choice has b een the de- 
fault for some time because calculation^SHIIZl have shown 
that revPBEapS! does not give rise to spurious binding. 
This choice of revFEE^; is thus conservative in that it 
ensures that all the improvement in binding compared to 
GGA calculations comes from the improvement in corre- 
lation ([I]). 

However, it has also been shown that revPBEa, is overly 
repulsive in the binding region of many vdW-bonded 
systems.!^ In the present paper our aim is to study the 
binding of the layers of V2O5 by qualitatively examining 
the improvement in binding in the vdW-DF calculations 
over the usual GGA calculations. We do not here aim for 
a "perfect" description of the binding in the sense of best 
possible numerical fit to experiment. Rather, by showing 
the results of using a few choices of sound exchange fla- 
vors, we illustrate the sensitivity to those choices and find 
the range of values of c and Eb within which we expect 
the physical values to be. A forthcoming publicatiorP^ 
addresses and discusses choices of exchange flavors with 
vdW-DF in extended (layered) systems. In the present 
presentation of V2O5 binding we pick two promising can- 
didates for exchange to use with vdW-DF, besides the 
well known PBEj; and revPBE^; functionals. These are 
C09^ (Ref. 41) and PW86:r (Ref. 42). 

GGA exchange may be described by the enhancement 
over LDA exchange in terms of the enhancement factor 

E.,^ j dvn[v)e'f'^[n{v))F^{s{v)) (5) 

where e^°^ = -c^n^/^ with = ?,^/^ / [An^'^) is the 
LDA exchange energy density and s = Cs|Vn|/n'*/'^ with 
Cs = l/(27r^/'^3^/'^) is the reduced electron density gradi- 
ent. The Fx{s) curves for the exchange flavors considered 
here are shown in Fig. [2] for a range of s- values. 

Dense materials have significant electron densities in 
most regions of the material and relevant values of s for 
those systems lie mainly in the rang^^ < s < 3. For 
sparse materials, with regions of very low electron den- 
sities, values of s > 3 also play a role and cannot be 
ignored. As is clear from Fig. [2] the various fiavors of 
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FIG. 2: (Color) Exchange functional enhancement factors fa; 
as functions of the reduced density gradient s ~ |Vn|/n^''^. 
Shown are the exchange flavors used here. 



exchange give different contributions at large values of s 
and the choice of exchange flavor is therefore more im- 
portant for sparse than for dense matter. 

In practice we include a particular flavor of exchange 
functional into the vdW-DF by subtracting the PBE^r 
energy contribution E^^^ from the total energy E'^^^ 
of the sc-GGA calculations (because the sc-GGA calcu- 
lations are carried out with the PBE flavor of GGA) and 
adding the value of the relevant exchange energy, calcu- 
lated from the sc-GGA electron density n. 



IV. RESULTS AND DISCUSSION 

A number of previous DFT studies of the V2O5 bulk 
system have used GGA calculations for obtaining the lat- 
tice constants of the orthorhombic unit cell.'Sitl2 The in- 
plane lattice parameters, in the present paper called a 
and h (but in parts of the literature called a and c), are 
found to be close to the experimental values. Table|l]lists 
for comparison the results of two experiments as well as 
the GGA results of Ref. HH The inter-plane binding dis- 
tance is the same as the third lattice parameter (here c, in 
parts of the literature 6) . For c a number of groups have 
found significantly larger values than the experimental 
value.-^— Ganduglia-Pirovano and Sauei^ carried out 
GGA calculations using PW91>S1 with PAW, emphasis- 
ing convergence and accuracy. This resulted in a c lattice 
constant even further away from the experimental value 
than the values from previous GGA calculations. 

The deficiency in the description of the binding from 
GGA calculations has been attributed to the presence 
of vdW interactions between the V2O5 layers.fiS] It is 
well known that GGA does not adequately describe vdW 
interactions,'^ thus if the layers of V2O5 are mostly 
bound to each other by vdW forces it is not surprising 
that GGA calculations do not give physically reasonable 



results. In the following we report our results from GGA 
calculations, from vdW-DF calculations, and discuss the 
effect of exchange and the binding character of V2O5. 



A. GGA results 

Exploring the V2O5 bulk system first with the 
(oACAPO-based) sc-GGA calculations we find, similar to 
the studies mentioned above, values of the a and h lat- 
tice constants in good agreement with experiment and a 
value for c about 11% too large when compared to exper- 
iment (Table |l]). The PBE total energy curve is shallow 
in the c direction. In order to determine the minimum 
with any reasonable accuracy we first find the approx- 
imate minimum, and then evaluate the total energy at 
125 points in (a, 6, c) parameter space in a small region 
around the approximate minimum. A three-dimensional 
polynomial fit to those points yields an energy minimum 
for the values of a, 6, and c as listed in Table|T]along with 
the results of our supplemental GPAW calculations. 

The shallow minimum of the PBE total energy curve 
at c = 4.87 A has the value 0.18 eV per unit cell com- 
pared to the layers taken apart, c — ?> 00 (c = 4.89 A with 
0.21 eV per unit cell for our GPAW calculations). This 
is in good agreement with the PW91 results of Ref. [12] 
where c — 4.84 A and a similar small surface energy in 
the cleavage plane was found. For truncated bulk they 
found the surface energy 0.048 J/m^, approximately cor- 
responding to a binding energy per unit cell of 0.25 eV. 
Our GGA results are also in good agreement with the 
PAW-based PBE results of Ref.HH In all four GGA cal- 
culations, the value of the binding energy is unphysically 
small, illustrating the deficiency of GGA. 



B. vdW-DF results 

With vdW-DF we find total-energy curves (at fixed 
values of a and 6) with enhanced binding compared to 
GGA (Fig. [3]). The position of the minimum of the total- 
energy curve depends on our choice of GGA exchange 
flavor, as discussed in the previous section. The optimal 
values of c and the binding energies for each choice of 
exchange functional are listed in Table [ij 

For vdW-DF with 009^; the value of c is close to the ex- 
perimental value, in fact the optimal value of c with C09 
is 2% smaller than the experimental value. PBEj., and 
PW862: as recommended in Ref. 1461 improve the binding 
over the results of GGA by decreasing c about 0.5 A, re- 
sulting in a value of c that deviates from experiment by 
only 2%. The binding energies range from 0.8 eV/unit 
cell (for revPBE^) to 1.3 eV/unit cell (PW862,). As ex- 
pected, revPBEa; leads to longer binding distances (9% 
larger than experiment) and a smaller binding energy 
than use of the other exchange flavors because revPBEj, 
is overly repulsive. 
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TABLE I: Equilibrium lattice parameters (a, b, and c) and 
interlayer binding energy per unit cell (Eb). Numbers in 
parenthesis are obtained by PW91 calculations and used in 
the calculations indicated. For the vdW-DF calculations the 
flavor of the exchange functional used is shown explicitly, 
with subscript x denoting the exchange part of the exchange- 
correlation functional. 







I) 










fAl 
[^1 


fAl 
L^J 




[pVl 


[vjrjr aj 


± 1 LLo WU 1 nj 












PBE with USPP 


11.52 


3.57 


4.87 


0.18 


24 


PEE with PAW^ 


(11.55) 


(3.58) 


4.89 


0.21 


21 


vdW-DF, revPBEj; 


(11.55) 


(3.58) 


4.72 


0.84 


23 


vdW-DF, revPBEo; 


(11.55) 


(3.58) 


4.72 


0.86 


22 


vdW-DF, PBE^ 


(11.55) 


(3.58) 


4.46 


1.28 


57 


vdW-DF, PW86a, 


(11.55) 


(3.58) 


4.46 


1.19 


61 


vdW-DF, C09, 


(11.55) 


(3.58) 


4.28 


1.10 


58 


Comparison 












PW91 with PAW'^ 


11.55 


3.58 


4.84 


0.25 




Experiment'' 


11.512 


3.564 


4.368 






Experimonf^ 


11.508 


3.559 


4.367 







''Using the DFT code GPAW; 

''revPBEa; extracted from the dacapo code; 

'^Ref. 1121 Ef, is here estimated as the cleavage energy, 

calculated from the given surface energy of truncated 

surfaces: 0.048 J/m^; 

''Ref. US 

''Ref.ll Powder of V2O5 (~ 420 nm crystals). 



%. -0.6 





1 1 1 1 1 

^ GGA(PBE) 




\ 




with revPBEx 




with C09„ - 




Wlin rVVoDv 




With PBEj^ 


\ \ exper. 




1 





4 4.2 4.4 4.6 4.8 5 5.2 5.4 

d[Al 

FIG. 3: ^vdW-DF ygjjjg various flavors of exchange. For 
revPBEa; we plot both the result of using the exchange en- 
ergy provided directly by DACAPO (dashed line) and using our 
external implementation (solid lines). The two stars indicate 
the optimal points of our two PBE-based GGA calculations 
and the vertical line the experimentally found value of c. 



to expand or compress) in the direction perpendicular to 
the layers, compared to other directions, the value of C33 
is also a good estimate of the bulk modulus of the system. 
Based on the likely most realistic exchange choices (i.e., 
excluding here the overly repulsive revPBEa;) we find that 
V2O5 bulk with C33 ~ 60 GPa is a somewhat more stiff 
material than graphite, which has"'^ C33 = 37-41 GPa. 



In Table H] we list two results based on revPBEa;. Pre- 
viously, we have used the values of the revPBE^ energy 
that were provided from sc-GGA DACAPO calculations 
for input to the vdW-DF calculations. In this work we 
rely on our own external implementation for the exchange 
energy calculations. The implementation differences be- 
tween DACAPO and our external code are minor, and as 
shown in Fig.|3]for revPBEa;, the differences in the results 
are also minor, and can be ignored. 

While the discussion of best and physically most rea- 
sonable choice of exchange flavor to go with vdW-DF 
is still open, it is clear from our total-energy calcula- 
tions that the vdW forces do indeed play a prominent 
role in the binding of the V2O5 layers. Earlier, another 
group has examined the role of the vdW- forces in V2O5 
by means of the semi-empirical DFT-D approach (Ref. 
lisp , arriving at a similar conclusion. 

The vdW-DF total energy curves in Fig. [3] have small 
wiggles, in particular in the expansion region (c > 4.6 A). 
We find that the wiggles mainly arise from the exchange 
part of the energy, as also seen in, e.g., Refs. I35II71 The 
minimum of the iJ'vdW-DF Q^j-ves is found by a polynomial 
fit to the points closest to the minimum. From this fit we 
also extract the corresponding values of the elastic coef- 
ficient C33 (Table |l|. Because V2O5 is more soft (easier 



C. Numerical noise and exchange functionals 

Our DFT calculations are based on a description of the 
valence electron wavefunctions by USPP. The USPP are 
not normconserving and tend to give "noise" in regions of 
very low electron density, yielding some small (unphysi- 
cal) negative values of n. In order to handle these small 
negative values of n, DACAPO replaces on the FFT grid 
points all values n < nfloor.dac by rifloor.dac before calcu- 
lating the exchange energy. The floor is given in terms 
of the Bohr radius ao, ?^floor,dac = 10"^" |e|/ao ~ 
\e\/A^ being a small but positive minimum value of the 
density. 

Points in space in which this floor is applied become, 
via the formula s ^ \Vn\/n'^^'^ , points with extremely 
large (unphysical) values of s. For exchange choices like 
PBEa; and revPBEa; with a constant asymptote of Fx{s) 
the effect of this replacement is small: Very large values of 
s {'3> 100) contribute much less to the energy integral ([s]) 
than do moderate values of s (« 5 — 10), because Fx{s) is 
basically constant at large s (small n) and the integrand 
n'^^^Fx thus vanishes for small n. However, for more "ag- 
gressive" exchange functionals, in terms of growth of F^ 
for large values of s, this may cause problems in material 
systems with large regions of small or vanishing electron 
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density ("vacuum"), such as our reference calculationsPsl 
This is further discussed in a forthcoming pubhcation.^'^ 
Here, we simply state that the modification used in our 
external implementation is to instead remove contribu- 
tions to the integral of ([5]) that come from points in space 
with n < nfioor = 10^ jej/A"^, and that the procedure 
is not very sensitive to the particular value of rifioor- 

At the time when dense, bulk-like systems were the 
main systems examined with DFT, relevant values of s 
were considered to be small, in the approximate rang^^U 
< s < 3. In sparse matter this is not always the case: in 
our sparse matter systems we routinely work with (phys- 
ical) values of s which are several orders of magnitude 
larger than this. In our V2O5 reference calculations, con- 
sisting of a single layer of V2O5 surrounded by vacuum, 
almost 50% of the spatial grid points give rise to s-values 
larger than 12 both in our USPP-based electron density 
from the dacapo calculations and in the pseudo part of 
the electron density of our supplemental GPAW calcula- 
tions. In comparison, the bulk calculations at close-to- 
experimental separation (c — 4.3 A) show mostly small 
values of s. We find that s = 2.7 is the largest value 
in the GPAW calculation for bulk V2O5 while only 0.4% 
of spatial points have large values (s > 12) in the corre- 
sponding DACAPO calculation. 

As a measure of the effect of removing all points with 
n < nfloor (including points with negative values of n) 
we calculate an error measure rj. This error measure is 
defined as the (absolute value of the) contributions that 
would have come from the excluded points, weighted with 
the unit cell size V = abc, 

where the difference indicated by A is between contribu- 
tions from bulk and reference calculations. We find that 
this measure is approximately constant over the range of 
layer separations studied here, with values ranging from 
0.5 meV/A3 for C09^ to 1.0 meV/A^ for PW86^. The 
fact that r] is approximately constant with c fits well with 
our expectation that spacial points with n < nfloor ap- 
pear in the vacuum region of the unit cell in the reference 
calculations, a region that expands approximately as fast 
as V. 

We note that the numerical problems faced here are 
similar to those one of us has previously documented for 
a molecular dimer.l^ Then, as part of an early vdW-DF 
study, we investigated the energy difference between a 
far-apart benzene dimer in a unit cell of a certain (large) 
size and two isolated benzene molecules each in a unit 
cell of the same size as the original unit cell. Whereas 
physically this energy difference should vanish, a small 
but finite contribution SE^^i appeared.l^ This unphysical 
difference (5i?ref would not converge with, e.g., planewave 
energy cut off or other convergence parameters and it 
was found to originate mainly from the exchange part of 
the energy. In dimer and adsorption studies this problem 
can be overcome by simple error cancellation using a fixed 



size of the u nit cell and moving the material fragments 
around.l^^^ In bulk calculations the unit cell size must 
necessarily change, and no such full error cancellation is 
possible, leading to a usually small but nonvanishing rj. 



D. Binding of layers 

The results of our vdW-DF calculations, regardless 
of choice of exchange functional, indicate that the vdW 
forces are indeed important in the V2O5 system. We have 
further analysed the system by extracting the change in 
electron density, An(r), that arises when moving an iso- 
lated layer of V2O5 into the bulk structure (atom posi- 
tions kept fixed) . This short qualitative analysis is based 
on the n(r) extracted from the sc-GGA calculations, but 
we expect no major deviations from the result of using 
an electron density based on sc vdW-DF calculations. 

Fig. [4] shows that only a small amount of electron 
charge is moved when the V2O5 layers bind. Most of 
this change is on the vanadyl O atom and the V atom. A 
Bader analysis,!^ using the algorithm described in Ref. 
[52I and discussed in Ref. |53l confirms this picture of very 
little change of charge: only about 0.1 electrons move 
from V to vanadyl O, the charge on all other atoms is 
almost unchanged. This is in agreement with results of 
Ref. 54] 

We agree that the static electron density n{r) may not 
be the best tool to interpret the binding character, as sug- 
gested and discussed in Ref. i55i Nevertheless, we believe 
that qualitatively the combined interpretation from Fig- 
ure|4]and from the Bader analysis is sufficient to conclude 
that no short-range binding type (like ionic or covalent 
bonds) can explain the main part of the binding of the 
V2O5 layers. We find that inclusion of dispersive inter- 
action is essential to complete the description of V2O5 
cohesion. 



V. SUMMARY 

In this paper we have studied the role of the van der 
Waals bonding inside a layered oxide by m eans of the 
(by now) well tested vdW-DF met ho d.'^^'^ In particu- 
lar we have underlined the importance of accounting for 
the long range interactions that characterize this bond- 
ing in order to recover a structure of V2O5 close to that 
found by experiments. We have shown that the flavor of 
the exchange functional chosen for the vdW-DF calcula- 
tions can give noticeable differences in calculated binding 
distances and energies and the definition of the most ap- 
propriate one is still a matter of debate. Specifically, 
in this paper we have tested two promising functionals 
to be used with vdW-DF: 009^, and PWSe^^. We have 
also analysed effects of numerical noise that may arise in 
studying sparse matter within DFT. This is important 
because there are large regions with low density, regions 
that may cause errors in the evaluation of exchange. Fi- 



7 




FIG. 4: (Color) Change in electron density, An(r), as a result of assembling the layers of V2O5. Left: Isosurfaces of values 
0.02 |e|/A^ are shown, red indicates gain of electrons (loss of charge). The box indicates the region of the plot shown in the 
right panel. Right: The electron density difference is shown in a slice through the system, through the positions of the V atom 
and the corresponding vanadyl O atom of the next layer. Color scales are indicated in the right-hand bar, in units of |e|/A^. 
Figure created with vmd.'^ 



nally, we have identified how this error varies with the 
choice of exchange flavor. 
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